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ABSTRACT 

We present the results of a numerical investigation of current-driven instability in mag- 
netized jets. Utilizing the well-tested, relativistic magnetohydrodynamic code Athena, 
we construct an ensemble of local, co-moving plasma columns in which initial radial 
force balance is achieved through various combinations of magnetic, pressure, and 
rotational forces. We then examine the resulting flow morphologies and energetics to 
determine the degree to which these systems become disrupted, the amount of kinetic 
energy amplification attained, and the non-linear saturation behaviors. Our most sig- 
nificant finding is that the details of initial force balance have a pronounced effect on 
the resulting flow morphology. Models in which the initial magnetic field is force-free 
deform, but do not become disrupted. Systems that achieve initial equilibrium by bal- 
ancing pressure gradients and/or rotation against magnetic forces, however, tend to 
shred, mix, and develop turbulence. In all cases, the linear growth of current-driven 
instabilities is well-represented by analytic models. GDI-driven kinetic energy amplifi- 
cation is slower and saturates at a lower value in force-free models than in those that 
feature pressure gradients and/or rotation. In rotating columns, we find that mag- 
netized regions undergoing rotational shear are driven toward equipartition between 
kinetic and magnetic energies. We show that these results are applicable for a large va- 
riety of physical parameters, but we caution that algorithmic decisions (such as choice 
of Riemann solver) can affect the evolution of these systems more than physically 
motivated parameters. 

Key words: instabilities - magnetohydrodynamics (MHD) - methods: numerical - 
galaxies: jets - pulsars: general 



1 INTRODUCTION 

Collimated, magnetized flows are ubiquitous in astrophysics, 
but detailed studies of the stability and evolution of such 
systems have focused primarily on the shear-driven Kelvin- 
Helmholtz instability (KHI, see e.g., Hardee 2006 or Salvesen 
et al. 2012 for review), where magnetic forces are often rel- 
egated to a supporting role. In fact, much of the pioneering 
work on the intrinsically magnetized family of current-driven 
instabilities (GDI) was conducted by the plasma physics 
community with the goal of understanding confined plasma 
columns in the context of Tokamak-type fusion experiments. 
Nevertheless, there are many obvious astrophysical circum- 
stances where GDI have the potential to play a significant role. 

Jets emerging from accretion onto black holes, for exam- 
ple, are very likely to be magnetically launched (e.g., Lovelace 
1976, Blandford 1976, Blandford & Payne 1982) and subse- 



quently magnetically dominated near their points of origin. 
Such flows will generally be subject to GDI, but the ramifi- 
cations of instability in this context remain unclear, particu- 
larly when the jet is bounded by a strong shear layer. Even 
if a shear layer confines lateral jet expansion and disrup- 
tion caused by GDI, one could still imagine GDI-driven dis- 
sipation resulting in variability signatures, such as the very 
rapid TeV variability recently observed in blazars Mrk 501 and 
PKS 2155-304 (Aharonian et al. 2007; Albert et al. 2007). At 
greater distances from their origin, jets from active galactic 
nuclei (AGN) are also inferred observationally to transition to 
a kinetic-energy dominated flow (Sikora et al. 2005), a pro- 
cess which itself may very well involve GDI-driven dissipation. 
Even in such a weakly magnetized regime, however, GDI have 
the potential to grow and subsequently to disrupt otherwise 
laminar flows. 

This is not an issue exclusively relevant to ultra- 
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relativistic jets. In the context of filled supernova remnants 
(i.e., plerions), notably including the Crab Nebula, GDI have 
been identified as a possible means of breaking axisymme- 
try in the outflowing pulsar-driven wind (Begelman 1998). If 
present, this process could obviate commonly made geomet- 
ric assumptions that constrain the magnetization of the flow 
at the wind shock to be significantly less than unity (e.g., Rees 
& Gunn 1974). Known as the "a problem," this inference has 
long been at odds with the fact that the wind is likely to be 
magnetically dominated near the pulsar (Goldreich & Julian 
1970) while simultaneously unlikely to convert enough of its 
magnetic energy before reaching the shock (e.g., Begelman 
& Li 1994), but GDI-driven asymmetry may provide a solu- 
tion to this seeming paradox. Furthermore, it is possible that 
dissipation events facilitated by instabilities such as the GDI 
could produce the recently observed, rapid, gamma-ray flares 
observed in the Grab Nebula by AGILE (Tavani et al. 2011) 
and Fermi (Abdo et al. 2011). 

There has been much discussion over the years in both 
the plasma physics and astrophysics communities on the lin- 
ear development of the GDI under various physical condi- 
tions, with much of the pioneering work emerging in the 
latter half of the last century. Lundquist (1951), for exam- 
ple, identified a criterion describing the ratio of poloidal 
to toroidal field that would cause a magnetized cylinder to 
become unstable. Kruskal & Schwarzschild (1954) further 
demonstrated analytically that a plasma column threaded by 
a purely toroidal field (known in fusion literature as a "Z- 
pinch") is characteristically unstable. This work was later gen- 
erahzed by Tayler (1957) to multiple mode numbers. Addi- 
tionally, Kadomtsev (1966) identified a criterion for which 
the axisymmetric m=0 mode of GDI (commonly referred to as 
the "sausage instabiUty") is stable. Useful overviews of the ba- 
sic physics of GDI, including distinctions between the internal 
(fixed-boundary) and external (free-boundary) kink modes 
and discussions of cylindrical column stability, have appeared 
in a number of textbooks, including Bateman (1978) and 
Goedbloed & Poedts (2004), and a review article by Boozer 
(2004). 

In the context of astrophysics, stability analyses have 
been conducted in the non-relativistic Hmit by a large num- 
ber of groups (Gohn 1983; Pietrini & TorricelH-Giamponi 
1989; GorbelH & TorricelH-Giamponi 1990; Appl & Gamen- 
zind 1992; Appl 1996; Appl et al. 2000; Kersale et al. 2000; 
Bonanno & Urpin 2011) under various sets of assumptions. 
Research investigating stability in the relativistic regime has 
been undertaken primarily in the force-free limit (Istomin & 
Pariev 1994, 1996; Lyubarskii 1999; Tomimatsu et al. 2001; 
Narayan et al. 2009), with the exception of Begelman (1998), 
which did not make the force-free assumption. This latter 
work will form the primary point of comparison between an- 
alytic theory and those of our simulations that include a pres- 
sure gradient that initially balances the magnetic forces. 

Magnetohydrodynamic simulations of GDI have been 
conducted in a number of different contexts. In the non- 
relativistic regime, Lery et al. (2000), Baty & Keppens (2002), 
Nakamura & Meier (2004), Li et al. (2006), Nakamura et al. 
(2007), Nakamura et al. (2008), Moll et al. (2008), and 
Ivanovski & Bonanno (2009) have explored both local and 
global models of GDI evolution. Relativistic numerical inves- 
tigations of GDI, however, have thus far been conducted by a 
relatively small number of groups. In a recent series of papers. 



Mizuno et al. (2009) and Mizuno et al. (2011) have run and 
analyzed local simulations of GDI both generally and in the 
specific context of pulsar wind nebulae. These models, sim- 
ulated in a reference frame that is co -moving with the bulk 
jet flow, were designed to explore the evolution of GDI un- 
der the most idealized circumstances. They further extended 
this work to local models bounded by a shear layer to ex- 
plore interactions between GDI and the KHI (Mizuno et al. 
2011), providing an extension into the relativistic limit of an 
earlier investigation by Baty & Keppens (2002) . Additionally, 
McKinney & Blandford (2009) have conducted a set of global 
general relativistic simulations, finding that their jets develop 
only limited substructure resulting from the action of GDI, al- 
though Mignone et al. (2010) identify "jet wiggling" in their 
global special relativistic jets that they do attribute to GDI. 

Our goal in this work is to construct a large number 
of idealized jet configurations to determine various circum- 
stances under which GDI develop and to examine the result- 
ing flow morphologies and energetics. Adopting a similar ap- 
proach to that of Mizuno et al. (2009), we will exclusively 
consider local models in which the observer is co-moving 
with the bulk flow, far removed from any jet boundaries. Un- 
like previous work, however, we will explore several differ- 
ent arrangements of the initial forces, considering systems in 
which the magnetic fields are force-free and those in which 
the J X B force is balanced by pressure gradients and/or ro- 
tation. The results of these simulations will be compared both 
to analytic estimates of GDI linear growth and to the results 
of previously conducted numerical simulations. Additionally, 
we will determine whether special relativity plays an impor- 
tant role in the evolution of these systems even in a reference 
frame for which the bulk velocity has been transformed away. 
Section 2 outlines our numerical methods and the initial con- 
ditions of our various models of magnetized plasma columns. 
Section 3 describes the results of our simulations, in terms 
of the simulated column morphologies, energetics, and how 
these quantities depend on various physical and numerical 
inputs. Finally, we discuss the implications of our models and 
present our conclusions in Section 4. 



2 NUMERICAL DETAILS 

The calculations presented in this work were conducted using 
Athena, a second-order accurate Godunov flux-conservative 
code for solving the equations of magnetohydrodynamics 
(MHD)^. The basic algorithms implemented in Athena are 
described by Gardiner & Stone (2005, 2008) with further de- 
tails (implementation and multi-dimensional tests) given in 
Stone et al. (2008). Specifically, we utiHze the dimensionally 
unsplit MUSGL-Hancock integrator ("Van-Leer") described by 
Stone & Gardiner (2009) combined with the constrained 
transport (GT) method of Gardiner & Stone (2005, 2008) to 
maintain the divergence-less nature of the magnetic field in 
multi-dimensions, extended to relativistic MHD (RMHD) as 
described by Beckwith & Stone (2011). 

The RMHD module within the Athena code implements 
a variety of Riemann solvers and spatial reconstruction meth- 
ods (as described in Beckwith & Stone 2011). Beckwith & 

^ The Athena code and a repository of test problems are maintained 
online at https : / / trac . princet on . edu/ Athena/. 
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Stone (2011) demonstrate that the choice of Riemann solver 
can play an important role in determining the spectral accu- 
racy of a numerical scheme. Even if a given scheme is stable 
and consistent, so that solutions converge to a weak solution 
of the conservation law, the solution itself may not be unique. 
That is, algorithms with different spectral performance may 
converge to different weak solutions of the conservation law 
As a result, it is important to assess the impact of different 
algorithmic choices (such as Riemann solver and spatial re- 
construction method) on the outcome of simulations of CDI. 
The impact of these choices, along with that of numerical res- 
olution, on the development of CDI is examined in detail §2.2. 



2.1 Initial Conditions 

In this work, we study the development of CDI in a set of ide- 
alized systems in order that we may understand the physics 
behind the phenomena that develop. To do so, we conduct 
three basic categories of simulation for which initial force bal- 
ance is achieved in three different ways. In the Newtonian 
approximation and in the absence of an external force, hy- 
dromagnetic force balance of an ideal MHD fluid can be ex- 
pressed as 



Table 1. Summary of Jet Column Simulations 



= -(vV)v-^ + -(JxB) 
P P 



(1) 



where v is the velocity field, p is the gas pressure, p is the 
density, B is the magnetic field, and J = V x B is the current 
density. In constructing this force balance equation, we have 
taken units where G = M = c = 47r = 1. If we further assume 
that the velocity field is restricted to a rotational field of the 
form V = Vcf){r)4> and consider only radial force balance, this 
expression simplifies to 



0=^-i^ + i(JxB)., 
r P P 



(2) 



where the first term is the (exclusively non-negative) centrifu- 
gal acceleration, the second is the pressure gradient, and the 
third reflects magnetic pressure and/or tension. In our first 
category of simulation, which we label force-free (FF), each 
of these terms is independently set to zero, meaning that the 
system initially features uniform pressure, no rotation, and a 
non-zero magnetic field. In the second type (PB) of simula- 
tion, the columns are initialized such that the pressure gradi- 
ent and magnetic forces balance one another in the absence 
of rotation. In the third type (RPB), initial force balance is 
achieved through a combination of non-zero rotational, pres- 
sure, and magnetic forces. Table 1 provides a listing of all 
production-level simulations, and we will now explain the dif- 
ferent types of initial conditions employed in each. In all of 
these simulations, we assume an adiabatic equation of state 
with a constant 7 = 5/3, except as noted. 

The magnetic fields in our models, all of which can be 
represented in cylindrical coordinates as B = B^c^ + B^z, 
fall into three generic categories that we refer to in Table 1 as 
sinusoidal, Komissarov, and force-free. For the sinusoidal field 
configuration. 



0.5Bo(l - cos[27rr/rc]) r < rc 
r > Tc 



(3) 



ID 


FIELD 


/3 


T(rA) 


OTHER NOTES 


PB 


sinusoidal 


0.3 


20 




PB-/33 


sinusoidal 


3 


20 




PB-/330 


sinusoidal 


30 


11 




PB-kom 


Komissarov 


0.3 


20 




PB-bzl 


sinusoidal 


0.3 


20 


Bz = —B(f) 


PB-bzO.2 


sinusoidal 


0.3 


20 


Bz = -0.2B^ 


PB-gl.33 


sinusoidal 


0.3 


20 


7 = 4/3 


PB-newt 


sinusoidal 


0.3 


20 


Newt. 


PB-vrand 


sinusoidal 


0.3 


20 


Vkick randomized 


RPB 


sinusoidal 


0.3 


69 


Newt. 


RPB-slow 


sinusoidal 


0.3 


49 


Newt., decreased v^f, 


RPB-xrot 


sinusoidal 


0.3 


52 


Newt., extended i;^ 


FF 


force-free 


0.25 


600 




FF-hot 


force-free 


0.25 


600 


p ~ lOp, 7 = 4/3 


FF-v6 


force-free 


0.25 


600 


Vkick = O.OQVA 


FF-V25 


force-free 


0.25 


600 


Vkick = 0.25'UA 


FF-newt 


force-free 


0.25 


600 


Newt. 


FF-vrand 


force-free 


0.25 


600 


Vkick randomized 



Here, FIELD is magnetic field configuration, as described in Section 
2.1. (3 is minimum value of (3 = 2p/B'^ on the grid. T is simulation 
duration in units of ta = Tc/va, where va is the initial maximum 
Alfven speed on the grid and rc is the initial column radius. Newt, 
refers to Newtonian physics (as opposed to SR). 



coordinate, and rc is the characteristic column radius. Typi- 
cally, Bz = for sinusoidal fields, corresponding to a mag- 
netic pitch angle P = rBz/Bd) = 0. The sole exceptions are 
the PB-bzl and PB-bzO.2 models, where a poloidal field is 
present such that Bz{r) oc B(p{r), with the constants of pro- 
portionality provided for each model in Table 1. What we re- 
fer to as the Komissarov- type field (e.g., Komissarov 1999) is 
a purely toroidal field in which 



^ , . . Bo{r/rc) r <rc 

Bo{rc/r) r > rc 



(4) 



where Bq is a fiducial field strength, r is the radial cylindrical 



Note that our field model drops off more quickly outside of 
Tc than that of Komissarov (1999) in order that the field will 
be greatly reduced in magnitude as it approaches the compu- 
tational grid boundaries. Additionally, we point out that the 
Komissarov-type configuration possesses a formally discontin- 
uous magnetic field derivative at r = a feature that we 
avoid in our other field models (for a discussion of surface 
currents in initial conditions, see GourgouHatos et al. 2011). 
Finally, the force-free case is similar to that used in, e.g., Appl 
et al. (2000) and Mizuno et al. (2009), featuring a combina- 
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tion of toroidal and poloidal components: 



designed to excite the fastest-growing GDI modes: 



j,2 _^ j,2^ ' 



B,{r) = -Bo- 



j,2 _|_ ^2^ 



(5) 



This field configuration corresponds to a constant magnetic 
pitch of P = — rc. In the FF models, an envelope function is 
also applied in the outer 10% of the computational box to en- 
sure that the magnetic field smoothly asymptotes to zero be- 
fore intersecting the grid boundaries in directions perpendic- 
ular to the column axis. It is worth explicitly noting here that, 
despite the fact that we provide these expressions in cylindri- 
cal coordinates for simplicity the simulations themselves are 
conducted on a Cartesian grid. 

In all cases, the columns are initialized with a spatially 
uniform density profile. The initial pressure profiles for the 
FF models are also spatially uniform, chosen such that the 
sound speed is ^ 0.14c (except in the case of model FF-hot, in 
which the sound speed is ^ 0.57c). The PB and RPB models, 
however, feature pressure profiles that vary radially. Specifi- 
cally, the pressure in these two types of simulation is derived 
from the condition of initial force balance (Eq. 2), given as- 
sumed magnetic field and velocity profiles. Despite the fact 
that many of the simulations themselves incorporate special 
relativistic physics. Equation 2 is sufficiently accurate for ini- 
tial force balance, which we tested by running the various ini- 
tial configurations without applying a perturbation to make 
certain that the equilibrium did not evolve significantly. Ad- 
ditionally, the pressure is scaled so that the minimum ratio 
of gas-to-magnetic pressure (/3 = 2p/B'^) achieves the value 
specified in Table 1. All models except PB-/33 and PB-/330 fea- 
ture /3min < 1, reflecting the fact that portions of the com- 
putational grid are magnetically dominated. Additionally, the 
maximum Alfven speeds for the fiducial models range from 
0.1 (PB) to ^; A ^ 0.3 (FF), the latter of which begins to 
enter into the relativistically strong field regime. 

Finally, the velocity profiles that we employ consist of 
a combination of large-scale rotational flows (v^) and small 
appHed perturbations {5vr,Sv(p). In all but the RBP models, 
there is no large-scale rotation, so = 0. In the RBP models, 
the rotation profile is of the form 



'^^(/)(^) = 0.5^;o(l — cos[27rr/rv]), 



(6) 



which is identical in form to the sinusoidal magnetic field pro- 
files. This profile is motivated primarily by desiring rotation 
that is both slowly increasing from zero near the center of the 
column (so that it is numerically resolved) and zero also near 
computational grid boundaries perpendicular to the column 
axis. In the fiducial model RPB, vq — 0.2c and = Tc, so that 
the magnetic field and rotation field are spatially co-located. 
In model RPB -slow, the velocity is reduced such that vq — 0.1c 
(for Tv — Tc). Finally, in model RPB-xrot, the rotation profile 
is extended radially such that = 1.5rc (with vq — 0.2c), so 
that there is some rotation exterior to the magnetized region 
of the column. We should note that the RPB models are run 
exclusively using Newtonian physics. This is primarily for ease 
of establishing the initial force balance, although we intend to 
explore the potential effects of special relativistic large-scale 
rotation and alternative rotation profiles (see, e.g. Moll et al. 
2008) in future work. 

The perturbation applied is an |m| = 1 type perturbation 



5vr{r^ (/), z) = 5vq exp(— r/4)[sin(?7) cos((/)) + cos(?7) sin((/))] 
5V(f){r^ (/), z) = 5v{) exp(— r/4)[cos(?7) cos((/)) — sin(?7) sin((/))], 

(7) 

where 77 = 2tiz/ z^ize with Zsize representing the size of the 
domain in the direction parallel to the column axis. The value 
of ^^;o is chosen such that the perturbation magnitude is of 
order 1% of the initial maximum Alfven speed va in all but 
models FF-v6 and FF-v25 (see Table 1). 

All of our models are simulated on Cartesian meshes 
of uniform cubic zone size at all locations on the grid. The 
computational domain size for the FF models is XsizeAc = 
Vsize/rc = 2zsize/rc = 32, where the column axis is ini- 
tially oriented along the z direction. In the PB/RPB models, 
however, we increase the resolution to ensure that the pres- 
sure/rotation profiles are adequately resolved. This has the ef- 
fect of reducing the domain size for a fixed number of zones, 
so that for the PB/RPB models the computational grid spans 
Xsize/rc = Vsize/rc = 2zsize/rc = 6.4. These domains are 
spanned by 256 x 256 x 128 zones in all of the production-scale 
models Hsted in Table 1, and we discuss tests of this resolution 
in §2.2. In all models, standard Athena outflow boundaries 
are employed in the x and y directions while periodic bound- 
aries are employed in In practice, the simulations are halted 
before there are substantial interactions with grid boundaries 
orthogonal to the column axis. 

Times in the remainder of this work will be given in units 
of a representative Alfven crossing time ta, which we define 
as TA = Vc/vA, where va is the maximum initial Alfven speed 
on the grid and rc is the initial column radius. While the local 
Alfven speed will obviously be a function of both position and 
time, we emphasize that our definition of ta is only directly 
indicative of the minimum Alfven crossing time at t=0. 

2.2 Spectral Accuracy & Numerical Convergence 

As described in Beckwith & Stone (2011) and outlined above, 
choices regarding the Riemann solver and the order of ac- 
curacy of spatial reconstruction (along with numerical reso- 
lution) can have important consequences for the interpreta- 
tion of non-linear MHD phenomena observed in numerical 
simulations. In the case of CDI in relativistic jets, the major- 
ity of prior simulations (e.g. Mizuno et al. 2009; McKinney 
& Blandford 2009) have been conducted using approximate 
HLLE-like solvers, where the structure of the Riemann fan is 
diffused (Harten 1983). For these solvers, discontinuities that 
are exact solutions of the Riemann problem will evolve via 
numerical diffusion, even when the flow is at rest (Mignone 
et al. 2009). This diffusion plays an important role in deter- 
mining the overall spectral accuracy of the numerical scheme 
(Mignone et al. 2009; Beckwith & Stone 2011). This numer- 
ical diffusivity can be reduced (and the spectral accuracy of 
the scheme improved) by including the physics of contact and 
Alfven waves in the solution of the Riemann problem, e.g. 
as in the five-wave approximate Riemann solver (which we 
term 'HLLD') initially described by Miyoshi & Kusano (2005) 
in the Newtonian regime and later extended to special rel- 
ativity by Mignone et al. (2009). If the physics of the prob- 
lem under consideration depends on the propagation of ei- 
ther of these wave families (e.g. the contact wave in the case 
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Figure 1. A comparison of magnetic field isosurfaces between the 
HLLD and HLLE solvers for otherwise identical CDI simulations at 
t = 10 TA- The resolution used in these models is a factor of two 
lower than any of our production runs to emphasize how different 
Riemann solvers can produce strongly divergent outcomes. 



of counter-propagating shear flows, or Alfven waves in the 
case of the magnetorotational instability), then use of Rie- 
mann solvers where these waves are diffused is likely to re- 
sult in incomplete physical interpretations of simulation out- 
comes, even though the solution itself is mathematically cor- 
rect. The growth of CDI is driven by the propagation of (modi- 
fied) Alfven waves, so CDI-driven outcomes are likely affected 
by choices regarding the Riemann solver. This point is illus- 
trated in Figure 1, which shows a side-by-side comparison of 
low-resolution (128 x 128 x 64 zones) test versions of model PB 
as conducted with the HLLD and HLLE solvers using Athena 
computed using second-order spatial reconstruction (see be- 
low) . Clearly, the CDI-driven flow morphology at this resolu- 
tion is dramatically affected by the choice of solver, particu- 
larly in regions of high magnetization that are conspicuously 
absent for the model computed with the HLLE solver, a result 
consistent with those of Beckwith & Stone (2011). 

To illustrate that this is not a problem exclusively con- 
fined to low resolutions, we also present a comparison of en- 
ergy evolution for the full-resolution PB models described in 
Section 2.1 and hsted in Table 1. While we will defer an ex- 
tensive discussion of CDI energetics to Section 3, it will suf- 
fice here to note that instability manifests itself as an initially 
exponential growth of kinetic energy (the so-called "linear" 
growth phase) that occurs at the expense of magnetic energy 
until the instability achieves a non-linear saturated state. Fig- 
ures 2 and 3 illustrate sample dependences of the kinetic en- 
ergy evolution on solver. In Figure 2, we see that solutions 
computed with the HLLD and HLLE solvers (again using 2nd- 
order spatial reconstruction) exhibit nearly identical growth 
rates during the linear phase of CDI development in model 
PB. Near the peak of the linear CDI development, however, 
the two solutions begin to diverge slightly until the end of the 
simulations, at which point the final kinetic energy for case 
computed with the HLLD solver is approximately 1.5 times 
that of the HLLE case. In the FF models, the difference be- 
tween solutions computed with the two solvers is significantly 
smaller, only ^ 17% over a much larger number of crossing 
times. Again, these differences are entirely attributable to the 
choice of Riemann solver for configurations that were other- 
wise identical. To summarize, improving the spectral accuracy 
of the Riemann solver leads to dramatic changes in flow mor- 
phology and energetics; in particular the efficiency at which 




5 10 15 20 
Time (rj 

Figure 2. Kinetic energy evolution for various solvers and reconstruc- 
tion orders for the fiducial PB model. Each model is independently 
normalized so that the total initial energy (TE+BE+KE) is unity. The 
combination of the HLLE solver and 3rd-order (3p) reconstruction 
leads to a final kinetic energy that is approximately half that of the 
HLLD solver with 2nd-order (2p) reconstruction. This is a clear illus- 
tration that algorithmic differences can produce divergent behaviors 
that, in most cases, have stronger effects than even the variation of 
physical model parameters. 




100 200 300 400 500 
Time (rj 

Figure 3. Kinetic energy evolution for various solvers and reconstruc- 
tion orders for the fiducial FF model. Each model is independently 
normalized so that the total initial energy (TE-I-BE-I-KE) is unity. The 
final energies differ by at most ~ 17%, which is significant but much 
less pronounced than in the PB model. 

magnetic energy is converted into kinetic energy. Since the 
HLLD approximate Riemann solver contains a more physi- 
cally complete description of the breakdown of the Riemann 
problem at the cell interface for the RMHD equations than 
the HLLE approximate Riemann solver, we adopt the HLLD 
solver as our standard of comparison and utilize it for all of 
the simulations presented in §3. 

While the algorithms implemented in Athena are over- 
all second-order accurate in space and time, extensive test- 
ing has revealed that use of third-order spatial reconstruc- 
tion can lower truncation error and increase the accuracy of 
the solution (Stone et al. 2008) . One might expect therefore. 
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Figure 4. A resolution study of energetics for the PB models, compar- 
ing the fiducial model PB (lines) with a model in which the physical 
resolution is doubled (glyphs) . Each model is independently normal- 
ized so that the total initial energy (TE-I-BE-I-KE) is unity The mod- 
els exhibit good convergence, suggesting that the evolution of small- 
scale structure does not dominate the energetics of the system. 



that the divergence in flow morphology and energy evolution 
between solutions computed using the HLLE and HLLD ap- 
proximate Riemann solvers can be ameliorated by improving 
the accuracy of spatial reconstruction from second to third 
order for the HLLE Riemann solver. The data of Figures 2 
and 3 show that, in the PB models, the asymptotic value of 
the kinetic energy using 2nd-order reconstruction is a fac- 
tor of ^ 1.4 times that of 3rd-order reconstruction, while for 
the FF models, the ordering is reversed, with 3rd-order re- 
construction producing kinetic energies roughly 18% greater 
than 2nd-order. From this, we conclude that while the use 
of third-order reconstruction in combination with the HLLE 
solver does reduce differences in flow energetics between Rie- 
mann solvers in the FF case, the more complex dynamics in- 
troduced when (e.g.) gas pressure plays a significant role in 
force balance are best captured by the use of the more phys- 
ically complete HLLD approximate Riemann solver, even at 
lower spatial accuracy. This result is consistent with expecta- 
tions derived from extensive comparisons between high-order 
(e.g. WEN05) schemes and second-order Godunov meth- 
ods, where it was demonstrated that, for non-linear prob- 
lems, both high order schemes and second-order Godunov 
methods converge at first order (Greenough & Rider 2004). 
These authors conclude that, while higher order methods de- 
liver higher accuracy solutions at fixed computational cost for 
smooth solutions, for non-linear problems second-order Go- 
dunov methods deliver higher accuracy solutions than higher 
order schemes (again, at fixed computational cost). As a re- 
sult, we adopt second-order spatial reconstruction for the all 
of the simulations presented in §3. 

The simulations presented here were conducted in ideal 
relativistic MHD, and, as a result, none of the solutions in the 
nonlinear regime can be regarded as 'converged'. For conver- 
gence, a physical dissipation scale (provided by either, e.g., a 
Navier- Stokes viscosity or Ohmic resistivity) would have to be 
included in the problem (Beckwith & Stone 2011). Such cal- 
culations are extremely challenging in relativistic MHD (see 
e.g. Takamoto & Inoue 2011) and are well beyond the scope 



Figure 5. A resolution study of energetics for the RPB models, com- 
paring the fiducial model RPB (lines) with a model in which the phys- 
ical resolution is doubled (glyphs) . Each model is independently nor- 
malized so that the total initial energy (TE-I-BE-I-KE) is unity The 
models exhibit adequate convergence with maximum deviations of 
approximately 5% near the onset of saturation and asymptotic differ- 
ences in the non-linear regime of approximately 30 — 35%. 
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Figure 6. A resolution study of energetics for the FF models, compar- 
ing the fiducial model FF (lines) with a model in which the physical 
resolution is doubled (glyphs). Each model is independently normal- 
ized so that the total initial energy (TE-I-BE-I-KE) is unity The models 
exhibit good convergence during the linear growth phase, but the ki- 
netic energy evolves differently at the ~ 30% level well after GDI 
saturation. 



of the current work. In addition, it is likely that in relativis- 
tic, magnetized astrophysical jets, the dissipation scale would 
be separated by many orders of magnitude from the small- 
est scales accessible by current computational resources. In- 
stead, convergence of the solution must be judged during the 
linear growth phase of GDI. Using the choices of Riemann 
solver and order of spatial reconstruction outlined previously 
(i.e. the HLLD approximate Riemann solver and second or- 
der spatial reconstruction), we compute simulations of the 
PB, RPB and FF cases described in §2.1 at both the fiducial 
resolution (256 x 256 x 128 zones) and at twice this resolu- 
tion (using otherwise identical configurations). The data of 
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Figure 7. Magnetic field isosurfaces for the fiducial PB model at four different times (specified in units of ta). All variants of the PB model are 
eventually subject to CDI, becoming completely disrupted over time. 



Figure 4-6 demonstrate that the linear growth phase associ- 
ated with each class of initial condition (PB: t < 7ta\ RPB: 
t < 15ta; FF: t < 75ta) is well converged at the fiducial res- 
olution of 256 X 256 x 128 zones, which we utiHze for the re- 
mainder of this work. While the subsequent, non-linear, stage 
of the evolution cannot be used to judge convergence of the 
models, it is instructive to examine the differences that arise 
at higher resolution. Model PB exhibits non-linear evolution 
that differs by < 5% between the fiducial and higher reso- 
lutions; both quantitative and qualitative conclusions drawn 
from the fiducial resolution during the non-linear evolution 



will therefore remain unchanged at higher resolution for this 
model. Models RPB and FF exhibit non-linear evolutions that 
differ by ^ 35% between the two resolutions at late times. In 
both cases, the qualitative non-linear evolution is not strongly 
affected by resolution, however. For the RPB model in partic- 
ular, the ratio of the magnetic-to-kinetic energy remains qual- 
itatively similar (i.e. approximate equipartition) during the 
non-linear stage. As a result, we conclude that while quanti- 
tative conclusions drawn concerning the non-linear behaviors 
of models RPB and FF will be subject to resolution effects. 
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Figure 8. Magnetic field isosurfaces for the fiducial rotating RPB model at four different times (specified in units of ta). As is the case for the PB 
models, the RPB models eventually become disrupted, although their field strengths remain generally higher than those of the PB models. 



qualitative conclusions regarding the physics at work in these 
models will be resolution independent. 



3 RESULTS 

3.1 Morphology 

To begin, we will discuss the distinct morphological evolu- 
tions that result from the three different classes of initial force 
balance. Here, we will focus primarily on comparing the fidu- 
cial models (PB, RPB, and FF), deferring a detailed discussion 



of the other models to our analysis of the system energetics. 
Figures 7-9 show snapshots of the evolution of magnetic field 
strength in models PB, RPB, and FF. The fields are shown 
as isosurfaces, which is sufficient to illustrate the degree to 
which the field structures bend and/or become disrupted as 
the systems evolve. 

Figures 7 and 8 show that both the PB and RPB columns 
are strongly disrupted, eventually showing clear signs of tur- 
bulent evolution. Model PB, in particular, features field struc- 
tures that emerge on a large variety of spatial scales, ranging 
from the size of the system to smaller scales that are obvi- 
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Figure 9. Magnetic field isosurfaces for the fiducial FF model at four different times (specified in units of ta). Unlike the PB/RPB models, the 
FF models are merely deformed, rather than completely disrupted. Although the system continues to evolve slowly after t ~ 100 ta, the gross 
morphology does not change dramatically after this time. 



ously not directly relatable to the initial |m| = 1 scale per- 
turbation. Additionally, the PB column begins to become dis- 
arranged quite rapidly, in less than 10 ta, and achieves some 
degree of turbulence by 20 ta. Considering that the system 
has also expanded to become larger than the initial column 
radius used to characterize ta, this is rapid evolution indeed. 

The RPB model, too, becomes disrupted, with notable 
differences from the PB model. One dissimilarity is a result 
of the system's rotation. Because the perturbation is initially 
non-axisymmetric, rotation strongly affects the orientation of 



the initial deformation and subsequent evolution. We also find 
that the RPB model appears to deform less rapidly than model 
PB. One might naively have assumed that the RPB model 
would only evolve more rapidly than model PB, given that 
the two systems have similar initial field configurations and 
that rotation provides an additional source of free energy in 
the RPB system. In detail, however, the two models also fea- 
ture different initial pressure gradients to satisfy Equation 2. 
This can further affect the magnetic field strength since it is 
only the minimum value of /3 that is constrained, meaning 
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Figure 10. Relating the direction of column deformation to the initial perturbation for the PB model (left), the RPB model (center), and the FF 
model (right). Solid lines show the angular orientation of the initial perturbation as a function of z, where the glyphs depict the center-of-mass 
(COM) and the average location of the magnetic field (COB) at late times. In all cases, the magnetic field and center-of-mass are close to 180 
degrees out of phase, which is not surprising since the field expands into a region of initially uniform density. For the PB and FF models, the 
magnetic field is clearly correlated with the direction of the initial perturbation, as we would expect. The RPB model features an evolved state 
that is shifted in phase as a result of the system's rotation. 



that this value can be achieved at different locations in differ- 
ent models. That the RPB and PB models evolve over differ- 
ent timescales simply illustrates that the characteristic Alfven 
crossing time alone may not always be perfectly indicative of 
the GDI evolution timescale, particularly if the system is un- 
dergoing significant evolution. We further note that the final 
field structure in model RPB does not appear to have accessed 
as large a range of spatial scales as model PB. This is again 
likely to be a result of the slower evolution of this system. 

Finally, Figure 9 illustrates that the FF model is deformed 
but not completely disrupted. Specifically, it is clear that the 
column reacts to the initial perturbation, but it is equally ob- 
vious that no smaller scales develop appreciable power as the 
system evolves. Furthermore, there is very little morpholog- 
ical change in the system after t ^ 100 — 200 ta, especially 
where the strongest field is located. This is obviously a com- 
pletely different regime of behavior from what was seen in 
those models that achieved initial force balance through a 
combination of magnetic and pressure forces. Whereas model 
PB becomes largely disrupted by 20 ta, model FF continues 
to be dominated by the |m| = 1 large-scale deformation even 
after evolving for 600 ta and shows no sign of turbulence. 

To demonstrate that we are in fact seeing these systems 
respond to the applied perturbation, we show in Figure 10 
the relationship between the deformation of the column and 
the direction of the applied perturbation as a function of the 
position along the column axis. Specifically, we look at the 
angular position of both the center-of-mass (COM), defined 
in a given direction as T^piXi/T^pi, and the analogous quan- 
tity for the magnetic pressure (COB), defined as T^Bfxi/T^Bf. 
Here, the Newtonian expressions for the density, position, and 
the magnetic field are adequate to provide a qualitative sense 
of whether or not the perturbation and subsequent evolution 
are aligned. For models PB and FF, Figure 10 illustrates that 
the orientations of the COB (asterisks) and the initial pertur- 
bation (lines) are well-correlated at all heights within in the 
column. In contrast, the COM (plus signs) is out of phase by tt 
radians in these two models. Together, these facts suggest that 
the magnetic fields move in the direction of the initial pertur- 
bation while the mass moves in the opposite direction. This 
is not terribly surprising, given that a deformed set of mag- 
netic field isosurfaces will generally expand to bound a larger 



volume than initially cylindrical isosurfaces. The density in- 
side such a region is observed to decrease accordingly as a 
consequence of this expansion (and the magnetic flux freez- 
ing captured in ideal MHD), leading to the center of mass 
becoming displaced in the opposite direction. In practice, this 
process leads to only a small overall displacement of the cen- 
ter of mass because so much of the mass in the outer portions 
of the grid is located in regions of weak or zero magnetiza- 
tion and subsequently does not move. In model RPB, we see 
similar behavior as in model PB, but shifted out of phase as 
a result of the system's rotation. In this model, different time 
snapshots would naturally feature different phase shifts. 

3.2 Energetics 

To understand how and why these different morphologies 
evolve, we now examine the details of energy flow in these 
systems. Our primary goal is to determine the degree to which 
the magnetic energy is being converted into kinetic energy as 
would be indicative of CDI. Furthermore, we seek to deter- 
mine whether or not this behavior saturates and what levels 
of kinetic energy amplification are attained through the ac- 
tion of CDI. Before proceeding, we should define precisely 
what we mean when discussing the various forms of en- 
ergy. In practice, Athena treats the total energy, defined as 
phT'^ - p + |B|V2 + [|v|2|B|2 - (v . B)2]/2, where h is the 
relativistic enthalpy and F is the Lorentz factor (e.g., Mignone 
et al. 2007, Beckwith & Stone 2011). Subtracting off the rest 
mass energy, the kinetic energy is KE = (F — l)Fp. The mag- 
netic energy is taken to be all terms containing the magnetic 
field, such that BE = |B|V2 + [Iv^Bl^ - (v • B)^]/2. In 
our models, the velocity terms in the magnetic energy expres- 
sion are typically measured to be less than one percent of the 
total magnetic energy, indicative of the characteristically low 
velocities achieved in these flows. The thermal energy is the 
remainder after the rest mass, kinetic energy, and thermal en- 
ergies have been subtracted from the total energy. 

3.2.1 Fiducial Models 

Figures 11-13 show the energy evolution for models PB, RPB, 
and FF. In each case, the figures depict the average energy 
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Figure 11. Energy evolution over time for the fiducial PB model, nor- 
malized so that the total initial energy (TE+BE+KE) is unity, where 
TE= thermal energy, BE = magnetic energy, and KE= kinetic energy 
The exchange of magnetic into kinetic energy reflects the develop- 
ment of CDI, which appears to saturate after several ta- The cyan 
region indicates the approximate time range over which CDI experi- 
ence linear growth. 



Figure 12. Energy evolution over time for the fiducial RPB model, 
normalized so that the total initial energy (TE-I-BE-I-KE) is unity In 
this case, the total kinetic energy (which includes rotational energy) 
actually decreases and kinetic/magnetic energy follow the same trend 
after a few tens of ta, both decreasing as they dissipate into thermal 
energy 



contained in each form (kinetic, magnetic, and thermal) as 
measured over the entire computational grid and normalized 
such that the total initial energy is unity. Although some en- 
ergy exits through the open computational boundaries at late 
times as these models evolve, this normalization reflects the 
total energy measured at any given time with an accuracy 
of better than two percent for all of our fiducial simulations. 
For all three fiducial models, the initial energy distribution 
is dominated by thermal energy with the magnetic energy 
comprising at most a few percent of the total grid energy. Of 
course, one should keep in mind that this relative ordering is 
largely an artifact of the grid size convolved with the radial 
profile of initial magnetic energy. Clearly, magnetic energy can 
make up a substantially larger fraction of the total energy in 
local regions where the field is stronger, and in fact all fiducial 
models feature a minimum /3 < 1. In each model, the kinetic 
energy attributable to the applied perturbation is a very small 
fraction of the initial total energy, as can be seen in Figures 
11 and 13. The perturbation is equally small for model RPB 
(Figure 12), but that model also features a non-negUgible (at 
the 2% level) amount of initial kinetic energy reflecting the 
rotation profile. 

In Figure 11, we can see the basic pattern of energy ex- 
change that reflects the development of CDI. Specifically, the 
total kinetic energy on the grid in model PB increases by 
approximately four orders of magnitude over several Alfven 
crossing times. During that same time period, the magnetic 
energy begins to decrease. The kinetic energy increase begins 
to roll over at t 6.5 ta, achieving the peak amplitude at 
t ^ 7 TA. After this time, both the kinetic and magnetic en- 
ergies show a decrease that is attributable primarily to tur- 
bulent dissipation into thermal energy and some loss off of 
the computational boundaries orthogonal to the column axis. 
This picture is consistent with the morphological evolution of 
model PB (Figure 7) in which the column bends, increasing 
the total kinetic energy, and eventually reaches a turbulent 
state in which the kinetic energy increases are offset by dis- 
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Figure 13. Energy evolution over time for the fiducial FF model, 
normalized so that the total initial energy (TE-I-BE-I-KE) is unity. Al- 
though the CDI initially act to amplify the kinetic energy at the ex- 
pense of magnetic energy, the overall levels of kinetic energy are al- 
ways far below those of the magnetic and kinetic energy (unlike the 
fiducial PB and RPB models). The cyan region indicates the approxi- 
mate time range over which the CDI experience linear growth. 

sipation into thermal energy. The cyan region of Figure 11 
depicts the approximate range of time over which the CDI un- 
dergo linear growth. In this region, we measure the maximum 
growth rate of the instability to be approximately 2.0 r^^ for 
the fiducial PB system. 

We can compare this growth rate to analytic estimates 
that have been constructed for physically analogous systems. 
Begelman (1998), for example, provides a framework for es- 
timating the growth rate of the toroidal field CDI under very 
general circumstances. By inserting the initial conditions of 
our PB model into this framework, we can compute the ex- 
pected CDI growth rate on a zone-by-zone basis within our 
computational grid to find the maximum. Combining the ex- 
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Figure 14. Energetics as they depend upon the details of rotation in RPB models, independently normalized so that the total initial energy 
(TE+BE+KE) is unity In model RPB-slow (center), the rotational velocity is halved from the case of the fiducial model RPB (left). The kinetic and 
magnetic energies in each run converge to similar values, even though the kinetic/magnetic energies in models RPB-slow and RBP are initially 
both different and differently ordered. In model RPB-xrot (right), the rotation profile is extended radially beyond the magnetized region, which 
lessens this strong connection between the kinetic and magnetic energies. While all three RPB models succumb to GDI, the global energy profiles 
depend upon the spatial partitioning of energy in the system. 



pressions provided in equations 3.18 and 4.8 of Begelman 
(1998) and inserting our initial values, w^e find that the max- 
imum expected growth rate in model PB is ^ 1.9 r^^. This 
is within five percent of the value measured from the en- 
ergy evolution of our simulation shown in Figure 11. More- 
over, this demonstrates that maximum growth is actually be- 
ing achieved in a large enough region of the grid that this rate 
dominates the energetics of the system. 

The energy evolution in the fiducial rotating system RPB, 
on the other hand, is very different than that of model PB 
even though both columns deform and develop turbulence. 
As Figure 12 illustrates, the magnetic energy first experiences 
a period of growth and is never reduced to the levels seen in 
model PB. This is consistent with the visual impression from 
Figure 8 that there are characteristically stronger fields in that 
model than appear in model PB (Figure 7). Moreover, after 
approximately 14 r a, the kinetic and magnetic energies in 
model RPB follow very similar trends, decreasing to nearly 
identical values. As we will soon see, the degree to which 
these quantities are coupled will depend upon the details of 
the rotation profile. In practice, it is challenging to derive an 
accurate growth rate for GDI in rotating systems from the en- 
ergy profiles alone since the kinetic energy due to rotation is 
simultaneously evolving in time. For that reason, we limit our 
current discussion of growth rates to the PB and FF classes of 
simulation. 

Finally, the energy evolution of model FF (Figure 13) is 
grossly similar to that of model PB despite the very different 
column morphologies that the two models develop. The pri- 
mary differences in their energy evolution are that the kinetic 
energy in model FF experiences a growth of only approxi- 
mately two orders of magnitude before saturating and that 
the peak amplitude of the kinetic energy is achieved at a much 
later time, around 110 — 130 ta. Additionally, we note that 
the growth depicted in Figure 13 appears to have two distinct 
phases, featuring a roughly constant slope from t = — 35 ta, 
whereupon it steepens until saturation. Moreover, the level of 
kinetic energy always remains a small fraction of the total 
energy, reaching only a tenth of one percent at maximum. 
Thus, we see that GDI do act to increase the kinetic energy in 
both the FF and PB configurations, but also that the levels to 
which the kinetic energy is amplified depend strongly upon 



the model parameters. The cyan region in Figure 13 again 
marks the period of linear growth. In the case of the fiducial 
model FF, the initial growth rate is 0.04 r^^, although the 
maximum growth rate at later times (but prior to saturation) 
is approximately twice this value. 

In the case of an initially force-free field, we can com- 
pare this growth rate to that analytically estimated by Appl 
et al. (2000), for example. This work provides a very simple 
closed expression for the growth rate in their Table 1, which 
is cjmax ^ 0.13^; A for force-free fields of constant pitch angle 
(identical to those of model FF). This translates in our case 
to cjmax ^ 0.04, which is identical to the initial growth mea- 
sured in our simulations. Again, this is approximately a fac- 
tor of two less than the maximum growth that is eventually 
achieved before saturation of the GDI, but it may be that the 
column evolution changes the Alfven speed sufficiently in por- 
tions of the column that higher growth is achievable at later 
times. An alternative interpretation is that this second phase 
of linear growth reflects the development of a different mode 
number, but the presence of such a mode is not easily inferred 
from examination of the flow morphology. It is important to 
note that both the morphology and linear GDI growth rate of 
our FF model appear to be very similar to those reported in 
Mizuno et al. (2009) for a nearly identical field model, which 
they refer to as their a = 1 case. 



3.2.2 The Effects of Initial Conditions 

To explore how the evolution of GDI depends upon the phys- 
ical parameters in our models, we first examine how the de- 
tails of rotation contribute to the system energetics in the 
RPB models. Specifically, we seek to determine the origin of 
the increased magnetic energy seen in Figure 12 that seem- 
ingly runs contrary to the usual GDI development. Figure 14 
shows a comparison of the fiducial rotating model RPB (left) 
with models RPB-slow (center), in which is decreased, and 
RPB-xrot (right), in which the profile is extended to larger 
radii. Model RPB-slow also shows a strong convergence of 
kinetic and magnetic energies despite the fact that, unlike 
model RPB, the system was initialized with more magnetic 
than kinetic energy. This suggests that the two forms of en- 
ergy are very likely to co-evolve whenever the GDI become 
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Figure 15. The evolution of kinetic energy over time as a function 
of /3 for models PB (solid line), PB-/53 (dotted line), and PB-/330 
(dashed line). For each model, the initial kinetic energy is inde- 
pendently normalized to unity. Scaling the evolution by the model- 
dependent Alfven time (ta) illustrates that the linear growth rate of 
CDI is relatively model independent, although the magnitude of the 
total amplification and turnover time both increase with increasing 
/3. 



active, somewhat independent of w^hich energy form is dom- 
inant. Although w^e have not attempted to identify the exact 
details of this energy transfer, one can imagine a mechanism 
not unlike the magnetorotational instability (Balbus & Haw^- 
ley 1991) operating to convert rotational energy to magnetic 
energy once the column has deformed sufficiently that the 
toroidal field can experience radial shear. Of course, this pre- 
supposes that the two forms of energy are distributed simi- 
larly. The rightmost panel of Figure 14 illustrates that a rota- 
tion profile extended beyond the magnetized region leads to 
parallel, but well-separated energy evolution, as one would 
expect. 

We are also interested in determining how CDI develop- 
ment is affected by magnetic field strength as parameterized 
by the plasma [3. Figure 15 shows how the kinetic energy evo- 
lution in the PB series of models varies with three different 
field strengths, reflecting both subthermal and superthermal 
configurations. In all cases, the initial perturbation was set 
to the same (small) fraction of the initial maximum Alfven 
speed, meaning that the initial kinetic energies differed be- 
tween models. We have therefore independently normalized 
the initial kinetic energies to unity so that we may more easily 
cross-compare the models. Clearly, Figure 15 illustrates that 
the linear growth stages of the CDI are not dramatically af- 
fected by the initial field strength. The growth is not only par- 
allel but nearly identical for each model until approximately 
5 TA, where the models begin to differentiate from one an- 
other. Specifically, the maximum saturation level and satura- 
tion time both appear to increase with increasing /3 (i.e., re- 
duced field strengths) . 

Additionally, we would like to know the importance of 
the magnetic field profile in determining the evolution of 
the CDI. Figure 16 shows the evolution of the alternate 
Komissarov-type magnetic field that is counterbalanced by a 
gas pressure gradient. Both runs clearly show evidence of ki- 
netic energy enhancement as a result of CDI. As measured 
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Figure 16. Comparing the kinetic and magnetic energy evolution be- 
tween the fiducial PB model (lines) with model PB-kom (glyphs), 
which features a Komissarov-type field profile. Each model is inde- 
pendently normalized so that the total initial energy (TE-I-BE-I-KE) is 
unity. While both models are disrupted by the CDI, the details of the 
energetics, including the CDI growth rate, differ according to the ini- 
tial field geometry and energy distributions. The linear growth rates 
are comparable to those predicted analytically for the two configura- 
tions. 



from the kinetic energy evolution, the CDI growth rates in the 
two models are comparable but not identical, suggesting that 
the exact distributions of field and gas pressure do play a role 
in determining the system evolution, as the analytic theory 
presented in Begelman (1998) suggests. The growth rate in 
model PB-kom is measured to be ^ 1.39 r^^, which is about 
70% that of model PB. If we again model the expected growth 
rate based on the expressions provided in Begelman (1998), 
we anticipate a growth of cj ^ 1.27 r^^, which is within 10% 
of what is measured. The other significant difference intro- 
duced by our choice of field structure is that PB-kom achieves 
equipartition between magnetic and kinetic energies while 
the kinetic energy in model PB evolves to exceed equiparti- 
tion. We should avoid over-interpreting this feature, however, 
since the global equipartition shown in the figure will not nec- 
essarily be indicative of local energy distributions everywhere 
on the grid. Specifically, one can easily imagine that the differ- 
ent field distributions in the two models lead to different total 
magnetic energies even if both systems experience rather sim- 
ilar development of CDI. In fact, direct inspection of the data 
confirms that there are a few local regions in model PB-kom 
where the kinetic and magnetic energies differ by a few orders 
of magnitude even though this is not reflected in the global 
average. 

Considering that there are significant morphological and 
energetic differences between the FF and PB models, it is im- 
portant to examine the role of poloidal field in the evolution 
of these systems. Specifically, we seek to determine whether 
the presence of poloidal field is sufficient to retard the devel- 
opment of CDI when introduced into the PB class of models. 
Figure 17 shows the kinetic energy evolution for models PB- 
bzl and PB-bzO.2, in which \B^\ = l^^l and \B^\ = O.2|B0|, 
respectively. These models differ from all of the others in that 
the magnetic pitch of these field configurations is proportional 
to the radius, rather than being constant. First, we note that 
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Figure 17. Comparing the kinetic energy evolution between the fidu- 
cial PB model (line) with models featuring poloidal fields that are 
comparable to 100% (PB-bzl) and 20% (PB-bzO.2) of the toroidal 
field strength. Each model is independently normalized so that the 
initial kinetic energy is unity The linear growth of the GDI is not dra- 
matically affected, although the maximum kinetic energy achieved 
and saturation behaviors are dramatically different. 

the linear grow^th of GDI is quite similar in all cases, confirm- 
ing that the different grow^th rates seen betw^een the FF and 
PB models cannot be attributed to the presence of a strong 
poloidal field in model FF. That said, it is clear that the linear 
grow^th terminates early for PB-bzl, suggesting that, as ex- 
pected, poloidal fields do play a major role in the saturation of 
the GDI. Furthermore, the saturation levels of the instability 
decrease on average w^ith increasing poloidal field strength, 
although w^ith substantial fluctuations that are of order a few^ 
times the characteristic Alfven crossing time. Incidentally, the 
column morphologies that develop for the poloidal models 
are also intermediate between those of the strongly turbulent 
PB model and the less disordered FF class of simulation. 

Next, w^e explore the degree to v\Ahich the temperature 
of the column affects its evolution. Model FF-hot, in w^hich 
the gas pressure is a factor of ten greater than the rest mass 
energy, is unique among our models in that special relativity 
is essential for achieving a subluminal sound speed, even in 
the initial conditions. Figure 18 compares the evolution of 
kinetic and magnetic energy in FF-hot w^ith the fiducial FF 
model. Again, both models feature clear indications of the 
development and saturation of the GDI. In fact, the energy 
evolution in model FF-hot runs roughly parallel to that of the 
fiducial model, w^ith a displacement that reflects the different 
initial fractional distributions of energy. This suggests that the 
temperature of the gas, even as it enters into the properly 
relativistic regime, does not strongly affect the development 
of the GDI. 

Finally, w^e examine the influence of the applied perturba- 
tion. The perturbation used in the vast majority of our models 
is a large-scale non-axisymmetric kick designed to excite the 
|m| = 1 mode of the GDI. Generally, this kick has a strength 
of approximately one percent of the initial maximum Alfven 
speed so that the system is excited, but not initially pushed 
too far aw^ay from the equilibrium state. Still, it is w^orth deter- 
mining w^hether or not the exact magnitude of this kick affects 
the evolution of GDI in a significant manner. Figure 19 depicts 



Figure 18. Comparing the kinetic and magnetic energy evolution be- 
tween the fiducial FF model (lines) and model FF-hot (glyphs), in 
which the gas is relativistically hot (i.e., p ~ lOp). Each model is in- 
dependently normalized so that the total initial energy (TE-I-BE-I-KE) 
is unity The two models feature similar CDI growth rates and parallel 
long-term evolution. 

the energy evolution of initially force-free models in v\Ahich 
the kicks were 6% (FF-v6) and 25% (FF-v25) of the maximum 
Alfven speed in otherwise identical simulations. The most ob- 
vious difference between these models naturally appears in 
the overall kinetic energy levels, which reflect the higher ini- 
tial kinetic energies associated with larger kicks. Otherwise, 
the evolution of GDI, particularly in the non-linear phases, 
runs roughly parallel for each different model, suggesting that 
the precise value of the kick is not particularly important for 
the evolution of the system. 

One may also explicitly test to determine whether the 
|m| = 1 mode is really the GDI mode that we should be ex- 
ploring. Linear theory tells us that this is under many con- 
ditions the fastest growing mode (e.g., Begelman 1998), but 
it is important to check that this mode is naturally excited 
by a more generic perturbation than is employed in most 
of our models. Figure 20 compares the fiducial PB and FF 
models with their randomly perturbed analogues, PB-vrand 
(left) and FF-vrand (right). Here, the random perturbations 
are generated by using a random number generator to pro- 
duce a random (in strength and orientation) velocity in each 
computational zone such that the RMS velocity magnitude as 
measured over the entire grid is comparable to the maximum 
kick applied in the fiducial models. This is an admittedly sim- 
ple test that would not, for example, be expected to converge 
numerically for increasing resolution, given that the pertur- 
bations are applied at the grid scale. Still, if GDI develop on 
large scales naturally from this type of perturbation, we can 
rest assured that our simulations reflect a fairly robust set of 
generic outcomes. 

In Figure 20, we see the telltale sign of GDI develop- 
ment, namely a clear increase in the kinetic energy on the 
grid. There are some interesting points to note, however. First, 
the PB and PB-vrand models (right) achieve rather similar 
magnetic energy evolution, but different kinetic energy pro- 
gressions. Specifically, the kinetic energy in model PB-vrand 
demonstrates a slightly different growth rate and saturation 
level from that of the fiducial model PB. This is not entirely 
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Figure 19. Energy evolution in the FF models as it depends upon perturbation magnitude for kinetic energy (left) and magnetic energy (right). 
Each model is independently normalized so that the total initial energy (TE+BE+KE) is unity Larger perturbations (e.g., FF-v25) obviously lead 
to larger initial/final kinetic energies, but all perturbations produce roughly parallel evolutions at late times. 




Figure 20. Energy evolution in the PB (left) and FF (right) models, comparing whether or not the system is subjected to an organized \m\ = 1 
perturbation (lines) or a random perturbation (glyphs). Each model is independently normaUzed so that the total initial energy (TE+BE+KE) is 
unity. Both models are obviously subject to the growth of CDI from random perturbations, but the PB model is more affected by the details of the 
perturbation than is the FF model, which is simply displaced in time. 



unexpected since the velocity field has time to modify the 
pressure (and field) profiles, both of w^hich can affect the CDI 
grow^th rate in pressure-supported columns, before the CDI 
have finished grow^ing. Additionally, model PB-vrand achieves 
approximate equipartition betw^een magnetic and kinetic en- 
ergies at saturation, w^hich is different from the PB model for 
w^hich the kinetic energy eventually exceeds the magnetic en- 
ergy globally In the FF and FF-vrand models, on the other 
hand, w^e see that the kinetic energies evolve quite similarly, 
w^hich is consistent w^ith the fact that the analytic grow^th rate 
estimate for force-free columns includes only the characteris- 
tic Alfven speed and radius, neither of w^hich should change 
dramatically as a result of a random perturbation. The only 
significant difference is an offset in time, but this is only nat- 
ural given that the |m| = 1 mode w^as not specifically excited. 



3.2.3 Including Alternative Physics 

In addition to exploring the influence of our choice of model 
parameters, w^e also seek to determine the significance of the 
included physics. For example, although it is not obvious that 
these systems experience a significant amount of compression 
as they evolve, it is w^orth determining w^hether a change of 
the adiabatic index has any noticeable effects. Figure 21 com- 
pares the evolution of fiducial model PB w^ith that of PB-gl.33 
in w^hich the adiabatic index of 4/3 corresponding to a rela- 
tivistic fluid is employed. There are essentially no important 
differences betw^een these tw^o models beyond an approxi- 
mately 10% deviation in the amount by w^hich the kinetic 
energy is amplified after 20 Alfven crossing times. This is a 
minor enough discrepancy that w^e can safely conclude that 
the value of the adiabatic index does not play a major role 
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Figure 21. Comparing the effects of the value of the adiabatic in- 
dex between the fiducial model PB (lines) with 7 = 5/3 and model 
PB-gl.33 (glyphs), in which 7 = 4/3. Each model is independently 
normalized so that the total initial energy (TE+BE+KE) is unity The 
model evolution is fairly insensitive to this parameter, with maximum 
deviations of approximately 10% in the kinetic and magnetic energies 
by the end of the evolution. 



Figure 22. A comparison between the fiducial model PB (lines) and 
its Newtonian analog, model PB-newt (glyphs). Each model is inde- 
pendently normalized so that the total initial energy (TE-I-BE-I-KE) is 
unity. The two models evolve in a similar fashion, although deviations 
can be seen, particularly in the kinetic energy where the relativistic 
value is ultimately 50% greater than that of the Newtonian case. 



in the evolution of GDI. One could imagine also investigating 
the effects of a variable Synge-type (Synge 1957) equation of 
state, but the relatively small ranges of temperature (a factor 
of 10, typically) achieved in a given GDI simulation suggest 
that this would not have a major effect on these models. 

A more fundamental physical difference may be expected 
to result from the incorporation of special relativistic physics. 
While our models are always initialized with low fluid veloc- 
ities, the Alfven speeds in these magnetized columns can be 
quite large. Specifically, the fiducial models all feature maxi- 
mum VA > 0.1, with the FF model having initial values as high 
as VA = 0.32. This is sufficiently large that it is not obvious 
how significant a role special relativity will play in regulat- 
ing this speed as the flow evolves. Furthermore, it is unclear 
whether the fluid velocities that develop as a result of GDI will 
be high enough to be strongly affected by the natural speed 
limit imposed by special relativity. 

Figures 22 and 23 show the evolution of the PB-newt 
and FF-newt models, respectively. It is initially surprising that 
the FF models are actually less sensitive to the inclusion of 
special relativity than are the PB models, despite their nomi- 
nally larger initial maximum Alfven velocities. Direct exami- 
nation of the data reveals that this is likely attributable to the 
larger flow velocities achieved in the PB models. When the PB 
columns are disrupted, small regions can attain velocities of 
up to nearly 0.14c. This is three times the maximum veloc- 
ity reached in the FF models and sufficient to be affected by 
the presence/absence of relativity. As a result, the fiducial PB 
model features a final kinetic energy that is nearly 1.5 times 
that of model PB-newt, while the FF models differ by only 
a few percent. This demonstrates that relativity does play a 
role in the evolution of these models, but also that special rel- 
ativistic effects do not dominate the energetics, as might be 
expected in simulations featuring relativistic shear layers or 
rotation profiles. 




100 200 300 400 500 
Time (rj 

Figure 23. A comparison between the fiducial model FF (lines) and 
its Newtonian analog, model FF-newt (glyphs). Each model is inde- 
pendently normalized so that the total initial energy (TE-I-BE-I-KE) is 
unity. The FF models evolve more similarly than the PB models, with 
only a few percent differences between relativistic and Newtonian 
physics. 



4 DISCUSSION AND CONCLUSIONS 

We have reported on an ensemble of local co-moving jet sim- 
ulations designed to explore the role of current- driven insta- 
bilities in magnetized plasma columns. Our most significant 
discovery is that the details of initial force balance have a 
dramatic impact on the resulting column morphology. Despite 
the fact that the force-free, pressure-supported, and rotation- 
supported magnetized columns under consideration are all 
nominally unstable to GDI, they clearly develop structures 
that are very dissimilar from one another. Specifically, our 
force-free columns undergo a relatively modest degree of 
deformation that slows dramatically when the GDI achieve 
non-linear saturation. The pressure-balanced columns, on the 
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other hand, shred and develop turbulent structures that facili- 
tate continued mixing throughout the evolution of the system. 
While the inclusion of rotation in pressure-supported systems 
adjusts some of the details, these models too are characteris- 
tically unstable and undergo significant mixing resulting from 
the combination of CDI and radial shear. 

Despite these grossly different morphologies, the gen- 
eral pattern of CDI energy evolution is similar for all mod- 
els. All systems experience a linear growth phase of the CDI, 
which we have found matches analytic predictions with a 
high degree of accuracy. The force-free fields generally expe- 
rience much slower rates of CDI development and less kinetic 
energy amplification than their pressure-supported counter- 
parts. Perhaps the most interesting energy exchange occurs in 
the rotating systems, where mixing ensures that spatially co- 
located kinetic and magnetic energies will evolve to approx- 
imate equipartition, regardless of their initial ordering. Con- 
sidering that most astrophysical winds and jets should pos- 
sess some amount of angular momentum, with AGN jets be- 
ing particularly likely to feature very rapid rotation near their 
source, this largely unrecognized channel of energy exchange 
could potentially be very important. 

Moreover, we have shown that these behaviors are fairly 
generic outcomes for a wide range of physical and numerical 
parameters. The strength and/or structure of the magnetic 
field, for example, makes much less difference than whether 
the initial force balance was achieved through force-free 
fields, pressure gradients, or rotation. Introducing a poloidal 
field component into a pressure-supported configuration does 
not affect the linear growth rate of CDI, although it does even- 
tually adjust the non-linear saturation properties. The adia- 
batic index or neglect of special relativity makes only a minor 
difference, although we have yet to merge special relativity 
and rotation in a single model. One expects that special rel- 
ativity will play a far greater role in affecting flows bounded 
by a strong shear layer and/or rotating at relativistic speeds. 

Numerical considerations prove to be some of the most 
significant factors in determining the non-linear states of 
these systems. Although the choice of Riemann solver and 
variable reconstruction order do not strongly impact the lin- 
ear growth phase of CDI, the peak kinetic energy attained 
and saturation levels are affected by algorithmic choices to a 
degree beyond most of the physical parameters that we var- 
ied. This serves as a useful reminder that not all numerical 
approaches are equally applicable to all problems and, specif- 
ically, that the HLLE solver should be avoided for problems in 
which it is important to resolve contact and rotational discon- 
tinuities. Given that the effective resolution of HLLD is greater 
than that of HLLE (e.g., Mignone et al. 2009), this is likely to 
be even more of an issue for global jet simulations that typ- 
ically employ lower physical resolutions over larger domains 
than their local counterparts. 

The different methods of force balance correspond to 
distinct regions of astrophysical parameter space. Force-free 
field configurations in pulsar winds or astrophysical jets, for 
example, are most likely to occur nearest to points of ori- 
gin. Simple arguments about flux conservation suggest that 
toroidal field will eventually dominate over poloidal field in 
any geometrically expanding jet or asymmetric wind (Begel- 
man et al. 1984). Furthermore, particles are more likely to 
be heated by dissipation downstream of their origin, with 
baryons in particular retaining this heat more easily than elec- 



trons or positrons. This provides a convergence of evidence 
that pressure-supported models are more relevant at greater 
distances from the source. 

That AGN jets are observed to be collimated over many 
orders of magnitude in physical scale, however, suggests that 
the disruptive behaviors seen in our pressure-supported and 
rotating columns do not completely dissociate real, astro- 
physical systems. One obvious difference between our current 
models and astrophysical jets is the presence of a confining 
medium that forms a shear layer with the jet. If nothing else, 
such a layer would certainly interact with and modify any CDI 
modes that intersected it, as has been seen in simulations by 
Baty & Keppens (2002) and Mizuno et al. (2011), for exam- 
ple. Under some circumstances, the shear layer would also be 
expected to develop KHI modes that would have the potential 
both to disrupt the shear layer and to interact with the exist- 
ing CDI modes. It is also possible, as has been suggested in 
the weak-field limit by Baty & Keppens (2002), that the KHI 
and CDI modes conspire to stabiHze the jet against disruption. 
It may even be that the details of jet rotation and the shear 
layer have a pronounced effect on whether or not jets are 
stable against CDI, which may reconcile the results of global 
simulations conducted by McKinney & Blandford (2009) and 
Mignone et al. (2010), for example. We will explore the in- 
teractions between the models we have presented here and a 
shear layer in a future work. 
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